SEAK_Analysis
Matthew R. Dunkle
2/19/2018
Abstract
This document is for the analysis of data collected in Southeast Alaska to assess variation in food web structure in streams across a gradient of hydrologic and physio-chemical conditions. This is a workind document and will be backed up on github at the following location: https://github.com/mdunkle/SEAK-Analysis.git
Fish Data
First thing I want to do is figure out the total biomass of fish captured in each sampling event.
## # A tibble: 455 x 19
## # Groups: ReachID, Date [15]
## ReachID StreamID Date Species Length Weight Temp Wwidth
## <fctr> <fctr> <date> <fctr> <int> <dbl> <int> <dbl>
## 1 MT_Upper Montana Creek 2017-07-21 Coho 39 0.60 12 8.9
## 2 MT_Upper Montana Creek 2017-07-21 Coho 45 1.30 12 8.9
## 3 MT_Upper Montana Creek 2017-07-21 Coho 44 1.10 12 8.9
## 4 MT_Upper Montana Creek 2017-07-21 Coho 42 1.00 12 8.9
## 5 MT_Upper Montana Creek 2017-07-21 Coho 37 0.50 12 8.9
## 6 MT_Upper Montana Creek 2017-07-21 Coho 79 0.56 12 8.9
## 7 MT_Upper Montana Creek 2017-07-21 Coho 84 0.73 12 8.9
## 8 MT_Upper Montana Creek 2017-07-21 Coho 44 1.00 12 8.9
## 9 MT_Upper Montana Creek 2017-07-21 Coho 38 0.60 12 8.9
## 10 MT_Upper Montana Creek 2017-07-21 Coho 54 1.90 12 8.9
## # ... with 445 more rows, and 11 more variables: BFWidth <dbl>,
## # Depth <dbl>, StreamType <fctr>, SamplingEvent <fctr>, X <lgl>,
## # Species.1 <lgl>, CPUE <lgl>, StreamID.1 <lgl>, ReachID.1 <lgl>,
## # number_fish <int>, biomass <dbl>
Plotting the length/weight relationship for Coho, Cutthroat, and DVs at the MT/McginConf
Here I’m visualizing lengths by species and reach
Coho Data
Next, I’m visualizing just the Dolly Varden data and showing the distribution as boxplots.
Dolly Varden Data
Next, I’m visualizing just the Dolly Varden data and showing the distribution as a boxplots. 
Length Distribution by Species and Sampling Event
Here, I’m using density ridge plots to show the distribution of lengths of each species by the two sampling events. The second plot shows the same data, but bracketed by stream type (glacial, clearwater (cw), brownwater (bw), and a combination of two (e.g. (CW/BW)) or all three (combined)). These distinctions are somewhat arbitrary, as most systems in SEAK show a combination of all three stream types and no one system is purely one type. It may be more useful to represent each of these as bins of conditions along a gradient of turbidity, flow variability, or DOC concentration.
cohodv_data = SEAK %>% filter(Species == c("DollyVarden","Coho")) %>% group_by(StreamID, SamplingEvent, Species)
cohodv=ggplot(SEAK, aes(x=Length, y=StreamID, fill=Species))+
geom_point(aes(color=Species), alpha=.5)+
geom_density_ridges(alpha=.5,scale=1,size=.00001, rel_min_height=0)+facet_wrap(~SamplingEvent)+
theme_cowplot()+theme(strip.background = element_blank(), legend.title = element_blank())
cohodv
SEAK_MeanWeight_Summary = SEAK %>% group_by(ReachID, StreamID, SamplingEvent, Species) %>% dplyr::summarise(mean=mean(Weight), na.rm=T)
ggplot(SEAK_MeanWeight_Summary, aes(x=SamplingEvent, y = mean, color=Species, group=Species))+geom_point()+xlab(NULL)+ylab("Mean Weight")+
facet_grid(~StreamID, switch="both")+
theme(axis.text.x=element_blank(), axis.ticks.x = element_blank(),
strip.background = element_blank(), strip.placement = "outside", strip.text.x = element_text(angle=90, vjust=1),
panel.spacing.x=unit(0, "lines"))+
geom_path(aes(x=SamplingEvent, y = mean, group=Species), arrow = arrow(type="closed"), size=1, position = position_dodge(width=1), alpha=.5)## Warning: Removed 16 rows containing missing values (geom_point).## Warning: Removed 10 rows containing missing values (geom_path).## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
Catch Per Unit Effort Across the Streams
Here, I load the data for catch per unit effort stored in the file SEAK_CPUE.csv. The associated figure shows the average catch across sites within a stream system during July and September sampling with an arrow showing the trend in catch rate between the two events. Important to note, Peterson Creek had the highest catch rate (and likely density) of any site in July and almost no juvnile salmonids in September. This could be due to hypoxic conditions or disturbance caused by adult salmon spawning in high densities in this system.
SEAK_CPUE=read.csv("SEAK_CPUE.csv", head=T)
SEAK_CPUE_Summary = SEAK_CPUE %>% group_by(StreamID, SamplingEvent, Species) %>% dplyr::summarise(mean=mean(CPUE))
print(head(SEAK_CPUE_Summary))## # A tibble: 6 x 4
## # Groups: StreamID, SamplingEvent [2]
## StreamID SamplingEvent Species mean
## <fctr> <fctr> <fctr> <dbl>
## 1 Fish Creek July Chinook 1
## 2 Fish Creek July CoastRangeSculpin 3
## 3 Fish Creek July Coho 9
## 4 Fish Creek July Cutthroat 0
## 5 Fish Creek July DollyVarden 2
## 6 Fish Creek September Chinook 0catch_trends=ggplot(SEAK_CPUE_Summary, aes(x=SamplingEvent, y = mean, color=Species, group=SamplingEvent))+
#geom_point()+
xlab(NULL)+ylab("Mean CPUE")+
facet_grid(~StreamID, switch="both")+
theme(axis.text.x=element_blank(), axis.ticks.x = element_blank(),legend.title=element_blank(),
strip.background = element_blank(), strip.placement = "outside", strip.text.x = element_text(angle=90, vjust=1),
panel.spacing.x=unit(0, "lines"))+scale_y_log10()+
geom_path(aes(x=SamplingEvent, y = mean, group=Species), arrow = arrow(type="closed"), size=1, position = position_dodge(width=1), alpha=.5)
catch_trends
save_plot("2017_CatchTrends.pdf", catch_trends, base_width = 8, base_height = 6)Periphyton Data
In this section, I visualize and analyze the data from periphyton samples taken during the 2017 pilot study.
Auxilliary Data
UAS Data from Jason Fellman showing the variables measured in 2006-2007.
## `geom_smooth()` using method = 'loess'## Warning: Removed 7 rows containing non-finite values (stat_smooth).## `geom_smooth()` using method = 'loess'## Warning: Removed 6 rows containing non-finite values (stat_smooth).## `geom_smooth()` using method = 'loess'## Warning: Removed 15 rows containing non-finite values (stat_smooth).
Modeling Simlulations
This was a practice for extracting data from Stella simulations and plotting in ggplot.
Making Maps for the proposal
Here I’m doing some spatial work to make sampling location maps. For some reason, this code is no longer working and needs some attention.
r, echo=F, message=F, warning = F} #library(“tmap”) library(“tmaptools”) library(“sf”) library(“leaflet”) library(“rio”) library(“maptools”)
wbdhu14 = “WBDHU14.shp” wbdhu14geo=read_shape(file=wbdhu14, as.sf=T) qtm(wbdhu14geo)
wbdhu8 = “WBDHU8.shp” wbdhu8geo = read_shape(file=wbdhu8, as.sf=T) qtm(wbdhu8geo)
wbdhu12 = “WBDHU12.shp” wbdhu12geo = read_shape(file=wbdhu12, as.sf=T) qtm(wbdhu12geo, “NAME” )
wbdhu10=“WBDHU10.shp” wbdhu10geo=read_shape(file=wbdhu10, as.sf = T) qtm(wbdhu10geo)
nhdflowline=“NHDFlowline.shp” nhdflowlinegeo=read_shape(file=nhdflowline, as.sf=T) qtm(nhdflowlinegeo)
alaskacoastline=“alaska_63360_py.shp” alaskacoastlinegeo=read_shape(file=alaskacoastline) pcreek=“Peterson Creek” wbdhu14geo=cbind(pcreek, wbdhu14geo) seakpoints=“SEAK-Sampling-Locations.shp” seakpointsgeo=read_shape(file=seakpoints) juneauwater=“tl_2015_02110_areawater.shp” juneauwatergeo=read_shape(file=juneauwater) head(juneauwatergeo, 50) mendenhallglacier=subset(juneauwatergeo, FULLNAME==“Mendenhall Glacier”) juneauicefield=subset(juneauwatergeo, FULLNAME==“Juneau Icefield”) akglaciers = “01_rgi60_Alaska.shp” akglaciersgeo=read_shape(file=akglaciers)
tm=tm_shape(wbdhu12geo, ylim=c(58.35, 58.55), xlim=c(-134.8, -134.38))+ tm_fill(“white”)+tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Lynn Canal”,])+tm_fill(“gray20”)+ #tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Eagle River”,])+tm_fill(“white”, lwd=2)+tm_text(text=“NAME”)+tm_borders(“black”)+ # tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Herbert Glacier”,])+tm_fill(“white”, lwd=2)+tm_borders(“black”)+ #tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Mendenhall Glacier”,])+tm_fill(“white”, lwd=2)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Fritz Cove-Frontal Lynn Canal”,])+tm_fill(“gray20”, lwd=2)+tm_borders(“gray20”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Tee Creek-Frontal Lynn Canal”,])+tm_fill(“gray20”, lwd=2)+tm_borders(“gray20”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“19010301060612-Mendenhall Peninsula”,])+tm_fill(“gray20”, lwd=2)+ tm_shape(alaskacoastlinegeo)+tm_fill(“white”)+ tm_shape(juneauicefield)+tm_fill(“dodgerblue”)+#tm_text(“FULLNAME”)+ tm_shape(mendenhallglacier)+tm_fill(“dodgerblue”)+#tm_text(“FULLNAME”)+ tm_shape(akglaciersgeo)+tm_fill(col=“dodgerblue”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Mendenhall River”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”, xmod=6, ymod=4)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Steep Creek”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”, xmod=3, ymod=2)+tm_borders(“black”)+ tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Herbert River”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Peterson Lake-Peterson Creek”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“pcreek”, xmod=1, ymod=-1)+tm_borders(“black”)+ tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Montana Creek”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”, xmod=0,ymod=-6)+tm_borders(“black”)+ #tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Windfall Creek”,])+tm_fill(“white”,lwd=2)+tm_text(text=“NAME”)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“McGinnis Creek”,])+tm_fill(“gray”,lwd=2, alpha=.0)+tm_text(text=“NAME”, xmod=7.5, ymod=2)+tm_borders(“black”)+ #tm_grid(projection=“longlat”, labels.size=.5, alpha=.75)+ tm_shape(nhdflowlinegeo)+tm_lines(“black”, alpha=.5)+ tm_shape(seakpointsgeo)+tm_symbols(col=“black”)+ tm_scale_bar()+tm_style_classic() tm
save_tmap(tm, “AKFieldsites.png”, width=3000, height=3000)
# function to obtain US county shape get_US_county_2010_shape <- function() { dir <- tempdir() download.file(“http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_us_050_00_20m.zip”, destfile = file.path(dir, “gz_2010_us_050_00_20m.zip”)) unzip(file.path(dir, “gz_2010_us_050_00_20m.zip”), exdir = dir) read_shape(file.path(dir, “gz_2010_us_050_00_20m.shp”)) }
obtain US county shape
US <- get_US_county_2010_shape()
split shape
US_cont <- US[!(US$STATE %in% c(“02”,“15”,“72”)),]
US_AK <- US[US$STATE == “02”, ] US_HI <- US[US$STATE == “15”,]
create state boundaries
US_states <- unionSpatialPolygons(US_cont, IDs=US_cont$STATE) tmap_mode(“plot”)
m_AK <- tm_shape(US_AK, projection = 3338) + tm_polygons(border.col = “grey50”, border.alpha = .5, breaks = seq(10, 50, by = 5)) + tm_layout(“Alaska”, legend.show = FALSE, bg.color = NA, title.size = 0.8, frame = FALSE) print(m_AK) ```